Ad Hoc(二)

  view

社区检测

上面我们已经有了每条船之间的熟悉度的公式，也就是每条船已经建立起来了联系，接下来我们就要使用已有的熟悉度，对现有的船进行社区的识别。

社区检测算法LOUVain

参考论文：《Fast unfolding of communities in large networks》
参考链接：算法原理
开发环境：anaconda
官方文档：API
语言：python
安装依赖包：conda install python-louvain

安装完成后，来运行一个实例

包名是community，但在pypi上引用python louvain

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import community as community_louvain
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import networkx as nx

# load the karate club graph
G = nx.karate_club_graph()

# compute the best partition
partition = community_louvain.best_partition(G)

# draw the graph
pos = nx.spring_layout(G)
# color the nodes according to their partition
cmap = cm.get_cmap('viridis', max(partition.values()) + 1)
nx.draw_networkx_nodes(G, pos, partition.keys(), node_size=40,
                       cmap=cmap, node_color=list(partition.values()))
nx.draw_networkx_edges(G, pos, alpha=0.5)
plt.show()

会显示图片：
在这里插入图片描述

1
community.best_partition(graph, partition=None, weight='weight', resolution=1.0, randomize=None, random_state=None)

使用Louvain启发式计算最大化模块化的图节点的划分（或尝试…）
这是最高模块化的划分，也就是由Louvain算法生成的树状图的最高划分。
返回 字典，分区(节点:对应的社区编号)，社区编号从0到社区数
系数

  • graph networkx.Graph
    分解后的networkx图
    对于neworkx.Graph图，这个我也是真的很无语了（哈哈哈）
    先来看一下官方的解读吧图相关的学习链接
    官方应该是提供了一部分的测试用例，来帮助你生成一部分图。
    但是我们要将自己的数据填充到里面去，我们就需要自己生成一个networkx.Graph图

    • 首先，创建一个空的图
      1
      G = nx.Graph()

    然后向空图中添加一些节点：

1
2
3
4
5
6
7
8
9
#单个节点的添加
G.add_node('1')
#通过一个节点集合添加
G.add_nodes_from(['2', '3'])
#添加没有权重的边
G.add_edge(n1, n2, object=x)
G.add_edges_from([('1', '2'), ('1', '3')])
#添加一个带有权重的边
G.add_weighted_edges_from([(1, 2, 0.125), (1, 3, 0.75), (2, 4, 1.2), (3, 4, 0.375)])
  • partition dict, optiona
    算法将开始使用这个节点分区。这是一个字典，其中键是它们的节点，值是社区
  • weight:str, optional
    图中用作权重的键。默认为“权重”
  • resolution:double, optional
    会改变社区的大小
    • randomize:boolean, optional
      将随机化节点评估顺序和社区评估顺序，以在每次调用时获得不同的分区
    • random_state:int, RandomState instance or None, optional (default=None)
      如果int，random_state是随机数生成器使用的种子；如果RandomState实例，random_state是随机数生成器；如果没有，则random number generator是使用的RandomState实例np.随机.

      将我们的数据分别填充进两个列表中去

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
shipid = open("/Users/gorge/Desktop/newshipData/out/r1_shipid.txt")
rela = open("/Users/gorge/Desktop/newshipData/out/r1_relation.txt");
links_data = []
for line in rela.readlines():
    dataline = line.strip('\n').split("[")[0].split(" ")
    famil=0
    if str(dataline[2])!='0':
        shijian = line.strip('\n').split("[")[1][:-1].split(",")
        #先算一个总的时间T
        T=0
        for a in shijian:
            if a!='':
              T=T+(int(a.strip(" ")))
        for x in shijian:
            if x!='':
                yy = int(a.strip(" "))
                famil=(famil+(yy**2/T))
        links_data.append((str(dataline[0]), str(dataline[1]), round(famil,2)))

nodes_data = []
G = nx.Graph()
dataline2 = shipid.readlines()[0].split(" ")
for id_q in dataline2:
    if str(id_q)!='':
        nodes_data.append(str(id_q))

将添加好的边以及边之间的关系进行添加

1
2
G.add_nodes_from(nodes_data)
G.add_weighted_edges_from(links_data)

使用算法进行社区的区分

1
partition = community_louvain.best_partition(G)

其他的一些参数设置

  • 补充一下字典中的初始化 函数
    1
    comm.fromkeys([str(i) for i in range(max(partition.values())+1)],[])

使用from keys 这个函数第一个参数为 默认的key值，第二个值为默认的value的值
在使用过程中发现，默认的value值如果为列表对象的话，初始值的设置为同一个对象，也就是说底层的指针，指向同一个地址

  • louvain中的模块度
    因为我们的实在需要一个评价社区划分的一个标准，经过查找资料发现，模块度可以用来评判社区检测的结果。
    因为louvain本来就是通过模块度的增益分社区的, 通过阅读community的源码，找到了对模块度的定义。
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70

    def modularity(partition, graph, weight='weight'):
        """Compute the modularity of a partition of a graph

        Parameters
        ----------
        partition : dict
           the partition of the nodes, i.e a dictionary where keys are their nodes
           and values the communities
        graph : networkx.Graph
           the networkx graph which is decomposed
        weight : str, optional
            the key in graph to use as weight. Default to 'weight'


        Returns
        -------
        modularity : float
           The modularity

        Raises
        ------
        KeyError
           If the partition is not a partition of all graph nodes
        ValueError
            If the graph has no link
        TypeError
            If graph is not a networkx.Graph

        References
        ----------
        .. 1. Newman, M.E.J. & Girvan, M. Finding and evaluating community
        structure in networks. Physical Review E 69, 26113(2004).

        Examples
        --------
        >>> import community as community_louvain
        >>> import networkx as nx
        >>> G = nx.erdos_renyi_graph(100, 0.01)
        >>> partition = community_louvain.best_partition(G)
        >>> modularity(partition, G)
        """
        if graph.is_directed():
            raise TypeError("Bad graph type, use only non directed graph")

        inc = dict([])
        deg = dict([])
        links = graph.size(weight=weight)
        if links == 0:
            raise ValueError("A graph without link has an undefined modularity")

        for node in graph:
            com = partition[node]
            deg[com] = deg.get(com, 0.) + graph.degree(node, weight=weight)
            for neighbor, datas in graph[node].items():
                edge_weight = datas.get(weight, 1)
                if partition[neighbor] == com:
                    if neighbor == node:
                        inc[com] = inc.get(com, 0.) + float(edge_weight)
                    else:
                        inc[com] = inc.get(com, 0.) + float(edge_weight) / 2.

        res = 0.
        for com in set(partition.values()):
            ss = (inc.get(com, 0.) / links) - \
                   (deg.get(com, 0.) / (2. * links)) ** 2
            print(ss)
            res += ss
        print(res)
        return res

我们可以通过将源码直接下载下来使用，为了保持文件的结构，我这里使用的是pycharm进行编辑的。
有了源码，其实我们也可以通过已经安装的community_louvain进行调用

1
2
cxc = community_louvain.modularity(partition, G)
print("社区的模块度为：",cxc)

这个函数一共有两个参数，上面源码中介绍的已经很清楚了
partition 为已经分出来的社区，

1
{'66660': 0, '55555': 1, '53995': 2, '12912': 3, '25417': 4, '68714': 5, '65773': 5, '50698': 6, '46344': 1, '29743': 1, '25510': 7, '11669': 3, '48048': 6, '59550': 8, '61259': 9, '46705': 0, '56825': 3, '14418': 3, '50140': 3, '66883': 3, '26371': 4, '47501': 1, '49980': 1, '31548': 6, ...}

partion是根据 {节点id : 社区名}进行标记社区的
G为输入的networkx graph图

1
2
3
4
5
6
7
G = nx.Graph()
其中G包括两部分组成：
节点的集合
['66660', '55555', '53995', '12912', ...]
节点之间关系的集合
[('66660', '31602', 15.0), ('66660', '66140', 1.95), ('66660', '41621', 10.14),...]
其中第一个元组代表的意思是 节点66660与节点31602的熟悉度为15.0

因此，在使用同一批进行

社区可视化

现在社区已经分出来了，使用python的库 networkx也可以画出社区关系图，但是画出来的图如下：
在这里插入图片描述
这个图是静态的，而且数据多了之后，可以发现区分的不是很明显，而且不能直观的表示出点的标签（这一部分我暂时没有探索出来，所以我就换了一个工具进行可视化处理）

Echarts

这是一个基于we b界面的可视化工具，在之前的记录中，也有提及过。这个画出来的模板是
在这里插入图片描述
可以看出来，这个图片看起来更友好一些，这里因为只能展示静态的图片，这个其实是一个动态的，当鼠标放在每个点上可以显示当前点的信息，而且更清晰一些。

工具

我这里使用的是webstrom 进行编辑的

需要引入的文件

1
2
3
<script src="./node_modules/js/echarts.min.js"></script>
    <script src="https://apps.bdimg.com/libs/jquery/2.1.4/jquery.min.js"></script>
    <script src="./node_modules/echarts/dist/extension/dataTool.js"></script>

  • 其中echarts.min.js的生成是由Echarts的官网进行定制的，我们可以登陆echarts的官网，然后点击下载按钮，选择在线定制，选择自己要定制的图表，然后在下面点击下载即可。
    在这里插入图片描述
    出现这个界面就表明定制成功了。

  • 对于query.min.js 大家可以使用在线的，直接使用上面的链接就可以了，这样就要保证在有网络的情况下使用，也可以去官网下载jquery.min.js的文件，这是j query所必须的依赖文件。

  • datatool.js 这个是我在下载echarts的包的时候就已经有了的。这个包的具体操作有点记不清楚了。很久之前做的。

数据格式

这个图使用的是json这个数据格式
解读一下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
{
"type": "force",
  "categories": [{
      "name": "HTMLElement",
      "keyword": {},
      "base": "HTMLElement"
    }]
   "nodes": [
    {
      "name": "AnalyserNode",
      "value": 1,
      "category": 4
    }]
    "links": [
    {
      "source": 169,
      "target": 137
    }
  ]
  }

我上面只是举了一个例子
可以发现先是categories 表示社区的类别，每个社区的名称是什么，name的作用是显示标签用的，base是用来点击社区的一个成员的时候显示用的。
nodes 是社区的点，name是社区成员的名字，value 是用来显示该成员的权重，category是用来显示该成员是属于哪一个社区的，该项的值从0 - n表示上面分类的index
links 是表示社区成员之间是不是有关系，source源，target 目标节点，同样的道理他们对应的值 为nodes的index的值。
因此，我们如果想要使用这个工具，我们还需要自己写一个数据的转换工具。

数据的转换工具

分类格式的转换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
comm = {}
#comm = comm.fromkeys([str(i) for i in range(max(partition.values())+1)],[])
for i in range(max(partition.values())+1):
    comm[i] = []
for key,values in  partition.items():
    comm[values].append(key)
# for key,values in comm.items():
#     print("第",key,"个社区,成员个数为",len(values),"成员为",values)

# 将数据转换为json格式
# 分类
categories = []
for i in range(len(comm)):
    sc = {}
    sc["name"]=i
    sc["keyword"]={}
    sc["base"] = i
    categories.append(sc)
print(categories)

节点格式的转换

1
2
3
4
5
6
7
8
9
# 节点的转换
json_nodes = []

for i in nodes_data:
    sc = {}
    sc["name"]= i;
    sc["value"]=1
    sc["category"]= partition[i]
    json_nodes.append(sc)

节点关系的转换

1
2
3
4
5
6
json_links = []
for i in links_data:
    sc = {}
    sc["source"] = ship_idmap[i[0]]
    sc["target"] = ship_idmap[i[1]]
    json_links.append(sc)

数据以json的格式进行存储

1
2
3
4
5
6
7
ww = open("/Users/gorge/Desktop/www.json","a")
jsonxx = {}
jsonxx["type"]="force";
jsonxx["categories"]=categories
jsonxx["nodes"] = json_nodes
jsonxx["links"] = json_links
ww.write(str(jsonxx))

数据可视化的时候的调整

  • 符号的调整
    经过上面的处理，现在我们已经有一个json格式的数据文件了，但是j son里面是不会出现单引号的，这时候我们的数据格式是：
1
'type': 'force', 'categories': [{'name': 0, 'keyword': {}, 'base': 0},

这个样子在json里面会报错的。我们在用代码写入的时候，python的字典、链表等都需要转换为字符串的格式进行写入，在这个转化的过程中，python会将所有的双引号以及单引号都转化为单引号。因此，我们写入的时候，是单引号。
这里我们就不用代码在转译了，我们就直接使用webStorm中自带的工具，Edit->Find->Replace 进行单引号与双引号的替换即可。

  • 颜色的处理
    使用模型自带的颜色是一共有9种颜色，在gexf文件中指定每个node的所属类别后，如果不为这些node自定义颜色的话，会在生成的画像中会为每一个类别默认设置一个颜色，每个类别的node拥有同一种颜色。
    可以看到明显有一个问题就是当你的分类数量大于9时，颜色会再次重复回去，这样对于分类任务来说是没法接受的。
    既然在生成关系图的设置中没有任何关于默认颜色的设置，那么这个颜色设定必然在自己所引用的echarts.min.js库中，所以只要耐心一点在echarts.min.js库中找到对应的部分，然后修改echarts.min.js库的源码这个问题就解决了。
    这个大约在16123行
    在网上找了一个颜色库，将默认的颜色库替换掉即可
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
color: [
                '#c23531', '#2f4554', '#61a0a8', '#d48265', '#91c7ae', '#749f83',
                '#ca8622', '#bda29a', '#6e7074', '#546570', '#c4ccd3', 'rgb(250, 235, 215, 1)',
                'rgb(240, 248, 255, 1)', 'rgb(0, 255, 255, 1)', 'rgb(245, 245, 220, 1)', 'rgb(255, 228, 196, 1)',
                'rgb(127, 255, 212, 1)', 'rgb(240, 255, 255, 1)', 'rgb(0, 0, 0, 1)', 'rgb(255, 235, 205, 1)',
                'rgb(0, 0, 255, 1)', 'rgb(138, 43, 226, 1)', 'rgb(165, 42, 42, 1)', 'rgb(222, 184, 135, 1)',
                'rgb(95, 158, 160, 1)', 'rgb(127, 255, 0, 1)', 'rgb(210, 105, 30, 1)', 'rgb(255, 127, 80, 1)',
                'rgb(100, 149, 237, 1)', 'rgb(255, 248, 220, 1)', 'rgb(220, 20, 60, 1)', 'rgb(0, 255, 255, 1)',
                'rgb(0, 0, 139, 1)', 'rgb(0, 139, 139, 1)', 'rgb(184, 134, 11, 1)', 'rgb(169, 169, 169, 1)',
                'rgb(0, 100, 0, 1)', 'rgb(169, 169, 169, 1)', 'rgb(189, 183, 107, 1)', 'rgb(139, 0, 139, 1)',
                'rgb(85, 107, 47, 1)', 'rgb(255, 140, 0, 1)', 'rgb(153, 50, 204, 1)', 'rgb(139, 0, 0, 1)',
                'rgb(233, 150, 122, 1)', 'rgb(143, 188, 143, 1)', 'rgb(72, 61, 139, 1)', 'rgb(47, 79, 79, 1)',
                'rgb(47, 79, 79, 1)', 'rgb(0, 206, 209, 1)', 'rgb(148, 0, 211, 1)', 'rgb(255, 20, 147, 1)',
                'rgb(0, 191, 255, 1)', 'rgb(105, 105, 105, 1)', 'rgb(105, 105, 105, 1)', 'rgb(30, 144, 255, 1)',
                'rgb(178, 34, 34, 1)', 'rgb(255, 250, 240, 1)', 'rgb(34, 139, 34, 1)', 'rgb(255, 0, 255, 1)',
                'rgb(220, 220, 220, 1)', 'rgb(248, 248, 255, 1)', 'rgb(255, 215, 0, 1)', 'rgb(218, 165, 32, 1)',
                'rgb(128, 128, 128, 1)', 'rgb(0, 128, 0, 1)', 'rgb(173, 255, 47, 1)', 'rgb(128, 128, 128, 1)',
                'rgb(240, 255, 240, 1)', 'rgb(255, 105, 180, 1)', 'rgb(205, 92, 92, 1)', 'rgb(75, 0, 130, 1)',
                'rgb(255, 255, 240, 1)', 'rgb(240, 230, 140, 1)', 'rgb(230, 230, 250, 1)', 'rgb(255, 240, 245, 1)',
                'rgb(124, 252, 0, 1)', 'rgb(255, 250, 205, 1)', 'rgb(173, 216, 230, 1)', 'rgb(240, 128, 128, 1)',
                'rgb(224, 255, 255, 1)', 'rgb(250, 250, 210, 1)', 'rgb(211, 211, 211, 1)', 'rgb(144, 238, 144, 1)',
                'rgb(211, 211, 211, 1)', 'rgb(255, 182, 193, 1)', 'rgb(255, 160, 122, 1)', 'rgb(32, 178, 170, 1)',
                'rgb(135, 206, 250, 1)', 'rgb(119, 136, 153, 1)', 'rgb(119, 136, 153, 1)', 'rgb(176, 196, 222, 1)',
                'rgb(255, 255, 224, 1)', 'rgb(0, 255, 0, 1)', 'rgb(50, 205, 50, 1)', 'rgb(250, 240, 230, 1)',
                'rgb(255, 0, 255, 1)', 'rgb(128, 0, 0, 1)', 'rgb(102, 205, 170, 1)', 'rgb(0, 0, 205, 1)',
                'rgb(186, 85, 211, 1)', 'rgb(147, 112, 219, 1)', 'rgb(60, 179, 113, 1)', 'rgb(123, 104, 238, 1)',
                'rgb(0, 250, 154, 1)', 'rgb(72, 209, 204, 1)', 'rgb(199, 21, 133, 1)', 'rgb(25, 25, 112, 1)',
                'rgb(245, 255, 250, 1)', 'rgb(255, 228, 225, 1)', 'rgb(255, 228, 181, 1)', 'rgb(255, 222, 173, 1)',
                'rgb(0, 0, 128, 1)', 'rgb(253, 245, 230, 1)', 'rgb(128, 128, 0, 1)', 'rgb(107, 142, 35, 1)',
                'rgb(255, 165, 0, 1)', 'rgb(255, 69, 0, 1)', 'rgb(218, 112, 214, 1)', 'rgb(238, 232, 170, 1)',
                'rgb(152, 251, 152, 1)', 'rgb(175, 238, 238, 1)', 'rgb(219, 112, 147, 1)', 'rgb(255, 239, 213, 1)',
                'rgb(255, 218, 185, 1)', 'rgb(205, 133, 63, 1)', 'rgb(255, 192, 203, 1)', 'rgb(221, 160, 221, 1)',
                'rgb(176, 224, 230, 1)', 'rgb(128, 0, 128, 1)', 'rgb(255, 0, 0, 1)', 'rgb(188, 143, 143, 1)',
                'rgb(65, 105, 225, 1)', 'rgb(139, 69, 19, 1)', 'rgb(250, 128, 114, 1)', 'rgb(244, 164, 96, 1)',
                'rgb(46, 139, 87, 1)', 'rgb(255, 245, 238, 1)', 'rgb(160, 82, 45, 1)', 'rgb(192, 192, 192, 1)',
                'rgb(135, 206, 235, 1)', 'rgb(106, 90, 205, 1)', 'rgb(112, 128, 144, 1)', 'rgb(112, 128, 144, 1)',
                'rgb(255, 250, 250, 1)',  'rgb(0, 255, 127, 1)', 'rgb(70, 130, 180, 1)', 'rgb(210, 180, 140, 1)',
                'rgb(0, 128, 128, 1)', 'rgb(216, 191, 216, 1)', 'rgb(255, 99, 71, 1)', 'rgb(64, 224, 208, 1)',
                'rgb(238, 130, 238, 1)', 'rgb(245, 222, 179, 1)', 'rgb(255, 255, 255, 1)', 'rgb(245, 245, 245, 1)',
                'rgb(255, 255, 0, 1)', 'rgb(154, 205, 50, 1)'
            ],

这样我们在去运行测试一下，就可以得到我们的社区分类的输出图形的可视化了
在这里插入图片描述
这样处理过后，我们的社区的可视化现在就相对很好了一些。

到此为止，对于社区检测的 关系图的输入，以及社区检测出来的社区的可视化分析我们就完成了，接下来我们需要做的是关系图的边的权重公式的优化，以及社区检测算法的优化。

测量社区分区算法分出来的社区的稳定性

现在我们一共有9天的数据（暂时还没有将更多的数据处理出来，因此我们就先用9天的数据分成3个3天进行测试社区的稳定性）

  • 社区的稳定性 定义为（同一批数据使用同一个社区划分算法，划分出来的社区 的个数相同，并且每一个社区的成员保持稳定）
    ｜这里的函数为 FileRead 中的div3d，控制函数为getFen3(String path,String targe) ｜

社区检测算法Infomap

算法简介
工具使用：在进行infomap工具安装的时候，网络上有很多的教程。

  • 第一种方式：因为我使用的是anaconda ，所以就先尝试了anaconda的安装方式安装连接，这个方式安装的时候，我没有安装成功。
  • 第二种方式：在线运行，这种方式在使用起来还是非常方便的，数据集可以直接在线上传，然后分区出来的结果可以在线可视化，缺点就是可修改性比较差
  • 第三种方式：使用pip安装，这种方式我还没有尝试。

    简单来介绍一下在线infomap的使用

    在这里插入图片描述
    点击进去，我们可以看到，这个界面，在使用的时候，只需要点击load network这个按钮将数据加载进去即可。
    数据格式可以是加权图，也可以是无权图

  • 带权图的格式
    第一行的数据表示为66660 节点与31620节点 存在边，边的权重为15.0

    1
    2
    3
    4
    5
    6
    7
    66660 31602 15.0
    66660 66140 1.95
    66660 41621 10.14
    55555 29743 6.0
    55555 49980 6.0
    55555 48640 363.0
    55555 42866 9.92

  • 无权图的格式
    第一行的数据表示为66660 节点与31620节点 存在边

    1
    2
    3
    4
    5
    6
    7
    66660 31602
    66660 66140
    66660 41621
    55555 29743
    55555 49980
    55555 48640
    55555 42866

接下来点击Run infomap这个按钮即可，下面显示框会出现运行报告

最左边的是open in networknavigator这个按钮，这个会在线显示分类的结果。
这个按钮下面有两个选项，可以选择clu 这个按钮，这个按钮是展示的最后最终的分类，还有一个按钮是Ftree这个按钮，这个会显示每一层的社区检测算法的结果。

最后的分类结果为
在这里插入图片描述
我们还可以点击每一个大的社区，查看这个大的社区下面小的社区
在这里插入图片描述
infomap检测的社区是分层次的。

  • Ftree 报告
    目前为止，对于这个社区检测工具的使用已经完成。但是我们现在还不清楚检测结果的Ftree报告
    这个报告主要分为两部分
    • 第一部分
      1
      2
      3
      4
      5
      # path flow name node_id
      1:1:1 0.0235485 "51031" 51031
      1:1:2 0.0213197 "25310" 25310
      1:1:3 0.020168 "26306" 26306
      1:1:4 0.0192412 "25771" 25771

第一行的意思是节点51031的节点id是51031，属于的社区层次为第一个社区中的第一号社区作为社区中的1号成员
现在flow的这个值，现在还没搞清楚

  • 第二部分为
    1
    2
    3
    4
    5
    *Links undirected
    #*Links path enterFlow exitFlow numEdges numChildren
    *Links root 0 0 24 7
    1 2 5.27517e-06
    1 4 6.39821e-05

enterFlow和exitFLow的值现在还没搞清楚

friendship 计算

通过大量阅读论文，找到了一篇跟我们场景比较像的论文《An Enhanced Community-based Routing Assisted by Ferry in Opportunistic
Networks》，我们对我们的friendship 的计算公式进行了改进。用两条船最后一次相遇的时刻减去第一次相遇的时刻，作为船的总的相遇间隔，然后乘二作为分子。分母是将每一次相遇的间隔的平方和累加。
然后对我们的计算公式，以及统计函数进行了改进。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import sys
ship = open("/Users/gorge/Desktop/newshipData/input/re50.txt")
info = open("/Users/gorge/Desktop/newshipData/r50port_ship.txt","a")
xx = ship.readlines()
T = 12957
res = []
count = 0;
for x in xx:
    v = x.strip("\n").split(" ")
    if len(v)>2:
        v = v[2:]
        SUM = 0
        if v[0]!='12960':
            #print(v)
            count+=1;
            for vi in v:
                SUM+=(((int(vi)**2)/2))
            if SUM ==0.0:
                SUM = T-3
            else:
                sx =T/SUM
                res.append(sx)
            info.write(str(sx)+"\n")
print(max(res))
print(min(res))
print(count)

绘制friendship的累积分布函数图

其中，最重要的是hist函数，其中familddd表示的是一个一维数组，也就是我们那要统计的信息，可以通过参数cucumulative来调节，默认为False，画出的是PDF，那么True画出的便是CDF直方图。PDF（figure1）可以观察到整个数据在横轴范围内的分布，CDF（figure2）则可以看出不同的数据分布间的差异性，也可以观察到整个数据的增长趋势和波动情况。
观察两种数据分布的差异，可能使用直方图就不是很直观，各种直方柱会相互重叠，我们只需更改直方图的图像类型，令histtype=‘step’，就会画出一条曲线来。
港口对应的船舶数目的分布关系图的统计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
ship = open("/Users/gorge/Desktop/newshipData/input/re50.txt")

xx = ship.readlines()
count=0;
port_shipc = {}
for i in range(47):
    port_shipc[str(i+1)] =0
print(port_shipc)
for x in xx:
    v = x.strip("\n").split(" ")
    if len(v)>3 and v[2]!='12960':
        port_shipc[v[0]] = (port_shipc[v[0]]+1)
print(port_shipc)
sd =0
for i in range(47):
    sd=sd+port_shipc[str(i+1)];
    tt = open("/Users/gorge/Desktop/newshipData/port_ship/"+str(i+1)+".txt","a")
    lii = "50 "+str(port_shipc[str(i+1)])+"\n";
    tt.write(lii)
print(sd)

下面使用matalab绘制的上述统计的关系图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
tic;
clear all;
data_dir = 'D:\MATLAB\port_ship\*.txt';
namelist = dir(data_dir)
len = length(namelist);
for i =1:len
    file_name{i} = namelist(i).name;
    x = load(file_name{i});
    H = plot(x(:,1),x(:,2));
    hold on;
end
legend("1号港口","2号港口","3号港口","4号港口","5号港口","6号港口","7号港口","8号港口","9号港口","10号港口","11号港口","12号港口","13号港口","14号港口","15号港口","16号港口","17号港口","18号港口","19号港口","20号港口","21号港口","22号港口","23号港口","24号港口","25号港口","26号港口","27号港口","28号港口","29号港口","30号港口","31号港口","32号港口","33号港口","34号港口","35号港口","36号港口","37号港口","38号港口","39号港口","40号港口","41号港口","42号港口","43号港口","44号港口","45号港口","46号港口");
LABEL_X = '通信半径/海里'; %X轴内容
LABEL_Y = '船舶数量/条';  %Y轴内容
%获取Y轴
hY = ylabel(LABEL_Y,'FontSize',15); 
hX = xlabel(LABEL_X,'FontSize',15);%同上
set(hY, 'FontSize',FONT_SIZE_LABEL);
set(hY,'Fontname','Times New Roman');
set(hY, 'Color', 'k');
set(hY, 'Interpreter','latex');
set(hX, 'FontSize', FONT_SIZE_LABEL);   %坐标轴名称大小
set(hX,'Fontname','Times New Roman');   %坐标轴字体
set(hX, 'Color', 'k');                  %坐标轴字体颜色
set(hX, 'Interpreter','latex');         %坐标轴格式

在这里插入图片描述

friendship的CDF关系图

1
2
3
4
5
6
7
tic;
clear all;
x=load('r50port_ship.txt')
[h,stats] = cdfplot(x)
title('Friendship Between Ships and Ports')
%xlim([0 0.0005])
legend('R50 CDF','Standard Normal CDF','Location','best')

在这里插入图片描述

船舶之间关系的统计

我们定义了一个船舶之间熟悉度的关系的公式，公式的含义是在一定的时间段下船舶相遇的平均等待时间，船舶的相遇定义为两个船之间的距离，小于一倍的通信半径。又在上面熟悉度的基础上进行了熟悉度的改进以及重新的统计计算。

船舶与港口之间关系的统计

船舶与港口之间的关系是将港口看作是一条船，然后计算每条船与港口之间的熟悉度进行定义的。

三个梯队船舶统计以及通信半径确定

  • 确定第一梯队 第二梯队 第三梯队的比例依次为3 ： 6 : 1
  • 第一梯队的船舶的确定方式为将港口看作是一条船舶的轨迹，然后计算船舶与各个港口的熟悉度，根据船舶回港口的次数， 统计回港口的船的数量以及船舶的ID 主要使用的文件为re5 表示的就是船舶与港口的关系
  • 第二梯队船舶的统计方式为计算出船舶与船舶之间的熟悉度关系，在所有的船号中，先将第一梯队的船号给记录下来，然后将第一梯队的每艘船统计相对应的相遇的船舶，然后在每一个第一梯队对应的相遇船舶中剔除掉，第一梯队的船舶。(如果不删除则要展现的是港口 - 对应第一梯队的船-对应第二梯队的船的相遇）
    具体的实现方法为：
    先使用python根据port_ship 的关系对应的文件是rex.txt ，统计出船与港口的关系，读取这一个文件，统计出第一梯队的船，保存到文件中，然后使用JAVA 语言读取第一梯队的文件，将船号保存到List中，并且使用Map使用第一梯队的船号作为Key ，Value设置为一个空的List ，进行初始化一个Map 集合，然后，读取船舶与船舶之间的关系的文件，rx_relation.txt 统计的船舶与船舶之间的关系，这个文件中包含三个信息，0 表示的是目标船A，1表示的目标船B，(1: ] 分别表示的是船A 与船B 等待下一次相遇的时间间隔。
    我们要按行遍历这个文件，首先要判断船A 与 船B 是不是全都属于第一梯队，这里使用and 判断，如果是船A 与船B中有一个是属于第一梯队 ，第二个不属于第一梯队，那么就找出不属于第一梯队的那个船，将该船号添加到Map的list中，这样就建立起了第一梯队的船与第二梯队的船的对应关系，这时候的第二梯队的船 是有重复船的。这样就可以将第二梯队与第一梯队的船建立起来了联系，写入到txt文件中。
  • 第三梯队的船 就是祛除掉第一梯队的船 与 第二梯队的船 剩下的船

    统计第一梯队的船舶的船舶id以及相应通信半径下的船舶的数量

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    ship = open("/Users/gorge/Desktop/newshipData/input/re50.txt")
    thepsn = open("/Users/gorge/Desktop/firstGroupShip/thepostship50number.txt", "a")
    thePs = open("/Users/gorge/Desktop/firstGroupShip/thepostship50.txt", "a")
    ships = open("/Users/gorge/Desktop/firstGroupShip/ships50.txt","a")
    xx = ship.readlines()
    count = 0;
    port_shipc = {}
    for i in range(47):
        port_shipc[str(i + 1)] = 0
    ship_portc = {}
    for i in range(47):
        ship_portc[str(i + 1)] = []
    shipc = {}
    for x in xx:
        v = x.strip("\n").split(" ")
        if len(v) > 3 and v[2] != '12960':
            port_shipc[v[0]] = (port_shipc[v[0]] + 1)
            ship_portc[v[0]].append(str(v[1]).split(".")[0])
            shipc[v[1].split(".")[0]] = 0
    for key, value in port_shipc.items():
        l = str(key) + " " + str(value) + "\n";
        print(l)
        thepsn.write(l)
    for key1, value1 in ship_portc.items():
        l1 = str(key1)
        for vx in value1:
            l1 = l1 + " " + vx;
        l1 = l1 + "\n";
        print(l1)
        thePs.write(l1)
    for key,vlaue in shipc.items():
        ships.write(str(key)+"\n")

数据可视化 Echarts 树状图json数据格式生成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
ship = open("/Users/gorge/Desktop/newshipData/input/re50.txt")
Treejson = open("/Users/gorge/Desktop/range50.json","a")
xx = ship.readlines()
count = 0;
port_shipc = {}
for i in range(47):
    port_shipc[str(i + 1)] = 0
ship_portc = {}
for i in range(47):
    ship_portc[str(i + 1)] = []

for x in xx:
    v = x.strip("\n").split(" ")
    if len(v) > 3 and v[2] != '12960':
        port_shipc[v[0]] = (port_shipc[v[0]] + 1)
        ship_portc[v[0]].append(str(v[1]).split(".")[0])

# 主要使用的是根据ship_portc进行拼接
# 首先创建一个父类 map
father = {}
father["children"] = []
father["name"] = "港口"
# 先将港口的名字拼接进去
for i in range(47):
    oil = {}
    oil["children"] = []
    oil["name"] = "港口"+str(i+1);
    for sx in ship_portc[str(i + 1)]:
        oil2 = {}
        oil2["children"] = []
        oil2["name"] = str(sx);
        oil["children"].append(oil2)
    father["children"].append(oil)
Treejson.write(str(father))

最后我们就得到了第一梯队与港口的关系,分别为5，10，20，50海里的关系。
在这里插入图片描述

统计第一梯队的船对应相遇的第二梯队的船

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
/**
   getTwoGroupShip("C:/Users/ouc/Desktop/firstGroupShip/ships5.txt","C:/Users/ouc/Desktop/result/r5_relation.txt","C:/Users/ouc/Desktop/twoGroupShip/one_twoGroupShips5.txt");

     * 
     * @param fistGroupPath 船舶与港口关系的文件 rex.txt
     * @param shiprePath 船舶与船舶之间关系的文件路径 rx_relation.txt
     * @param twog 第二梯队的船的保存路径
     */
	    public static void getTwoGroupShip(String fistGroupPath,String shiprePath,String twoGrouppath) {
	    	//首先是读取第一梯队的船 船号最为key list作为value
	    	FileRead tool = new FileRead();
	    	getTenMin wtool = new getTenMin();
	    	List<String> theFirstGroup = tool.readFileLine(fistGroupPath);
	    	//创建一个第一梯队对应的船的组
	    	Map<String, List<String>> fistmap = new HashMap<String, List<String>>();
	    	for (int i = 0; i < theFirstGroup.size(); i++) {
				fistmap.put(theFirstGroup.get(i),new LinkedList<String>());
			}
	    	List<List<String>> shiprela = tool.readFILE3(shiprePath);
	    	for(int i =0;i<shiprela.size();i++) {
	    		List<String> line = shiprela.get(i);
	    		if(Integer.parseInt(line.get(2))>0) {
	    			if(!line.get(2).equals("12960")) {
	    				String x1 = line.get(0).split("\\.")[0];
	    				String x2 = line.get(1).split("\\.")[0];
	    				
	    				//判断船A 与船B是不是都属于第一梯队
	    					if(theFirstGroup.contains(x1)&&!theFirstGroup.contains(x2)) {
	    						fistmap.get(x1).add(x2);
	    						System.out.println(x1+" "+x2);
	    					}
	    					if(!theFirstGroup.contains(x1)&&theFirstGroup.contains(x2)) {
	    						fistmap.get(x2).add(x1);
	    						System.out.println(x1+" "+x2);
	    					}
	    			}
	    		}
	    	}
	    	String value = wtool.WriteValue3(fistmap);
	    	wtool.writerFile(twoGrouppath, value);
	    }

统计第二梯队的船的数量以及ID

1
2
3
4
5
6
7
8
9
10
11
12
13
#统计第二梯队的船数目（重复的算一次）
ship = open("C:/Users/ouc/Desktop/twoGroupShip/one_twoGroupShips20.txt")
ship_id = open("C:/Users/ouc/Desktop/twoGroupShip/Twoships20.txt","a")
page = ship.readlines();
ff = {}
for x in page:
    if len(x)>0:
        xx = x.strip("\n").split(" ")[1:]
        for v in xx:
            ff[v]=0
print(len(ff))
for k in ff.keys():
    ship_id.write(str(k)+"\n")

统计第一梯队的船对应的第二梯队船的数目

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ship = open("C:/Users/ouc/Desktop/twoGroupShip/one_twoGroupShips20.txt")
ship_write = open("C:/Users/ouc/Desktop/twoGroupShip/one_twoShips20number.txt","a")
page = ship.readlines();
ff = {}
for x in page:
    if len(x)>0:
        cx = x.strip("\n").split(" ")[0]
        xx = x.strip("\n").split(" ")[1:]
        for v in xx:
            ff[cx]=len(xx)
for key,value in ff.items():
    s = str(key)+" "+str(value)+"\n"
    ship_write.write(s)
print(ff)
统计船舶与港口相遇次数 、每次相遇花费时间、以及完成相遇需要使用的时间的代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
//读取港口的文件，计算船与港口的伴随关系 
    public static List<String> getMeetportWith(List<String> firstline ,List<List<String>> second) {
		//定义一个存储上一个时刻的变量
		double pretime = 0.0;
		//定义一个变量用来 存储船的相遇次数
		int count = 0;
		String secshipid = second.get(0).get(1);
		String firshipid = firstline.get(0);
		System.out.println("港口"+firshipid+"相遇"+secshipid+"getMeetCount() start");
		List<String> lStrings = new LinkedList<String>();
		lStrings.add(firshipid);
		lStrings.add(secshipid);
		//初始化一个用来记录两条船时间间隔的链表
		List<String> banshui = new LinkedList<String>();
		//定义一个变量用来存储在船相遇的时候，伴随了几个间隔
		int b = 0;
		//定义一个旗帜来表示这段间隔是不是船开始伴随了
		boolean flag = false;
		//定一个变量来判断表示两个船是不是第一次相遇
		boolean firstmeetTime = true;
		//定义一个变量来存储第一次相遇的时刻
		double firtimeet = 0;
		//定义一个变量用来保存 最后一次相遇的时刻
		double endtimeet = 0;	
		if(second.size()>=4) {
		for(int i=0;i<second.size()-1;i++) {
			//获得第二条船的第i时刻的记录
			 List<String> secondline = second.get(i);
			 //记录一下当前时刻 
			 double curtime = Double.parseDouble(secondline.get(0));
			 //获得第一条船的第i时刻的 x y  
			 double firstx = Double.parseDouble(firstline.get(0));
			 double firsty = Double.parseDouble(firstline.get(1));
			 //获得第二条船的第i时刻的x y
			 double secondx = Double.parseDouble(secondline.get(2));
			 double secondy = Double.parseDouble(secondline.get(3));
			 //计算一下 两个点之间的距离 与 船舶的二倍的通信呢半径比较一下
			 double x = Math.pow((secondx-firstx), 2);
			 double y = Math.pow((secondy-firsty), 2);
			 double distance = Math.sqrt(x+y);
			 //设置一下船的通信范围
			 double contactRange = 5*1.87;
			 if(distance<=contactRange) {
				 if(firstmeetTime) {
					 firtimeet = curtime;
					 firstmeetTime = false;
				 }
				 endtimeet = curtime+3;
				 if(curtime-pretime>3) {
				 pretime = curtime;
				 count++;
				 flag= true;
				 b++;
				 }else {
					pretime = curtime;
					flag = true;
					b++;
					//System.out.println("b"+b);
				}
			 }else {
				 flag=false;
			 }
			 if(!flag && b!=0) {
				// System.out.println(flag);
				 banshui.add(String.valueOf(b*3));
				 b=0;
			 }
		}
		}
		if(b!=0) {
			//System.out.println("out"+b);
			banshui.add(String.valueOf(b*3));
		}
		//伴随的次数
		 lStrings.add(String.valueOf(count));
		 //每次伴随的时间
		 lStrings.add(String.valueOf(banshui));
		 //完成伴随需要花费的总的时间
		 lStrings.add(String.valueOf(endtimeet-firtimeet));
		System.out.println("船"+firshipid+"相遇"+secshipid+"次数为："+count);
		System.out.println("getMeetCount() end");
		return lStrings;
	}

统计船舶与港口相遇次数以及熟悉度相互对应的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import sys
ship = open("C:/Users/ouc/Desktop/result/re100.txt")
ship_number= open("C:/Users/ouc/Desktop/result/re100meet.txt");
info = open("C:/Users/ouc/Desktop/result/number_100rean.txt","a")
xx = ship.readlines()
yy = ship_number.readlines()
T = 12957

count = 0;
N_R = {}
for x in xx:
    v = x.strip("\n").split(" ")
    if len(v)>2:
        POR = v[0]
        shipId = v[1].split(".")[0]
        id = str(POR)+"_"+str(shipId)
        v = v[2:]
        SUM = 0
        if v[0]!='12960':
            N_R[id] = []
            #print(v)
            count+=1;
            for vi in v:
                SUM+=(((int(vi)**2)/2))
            if SUM ==0.0:
                SUM = T-3
            else:
                sx =T/SUM
                N_R[id].append(sx)
            #info.write(str(sx)+"\n")
for y in yy:
    vx = y.strip("\n").split(" ")
    if int(vx[2])>0:
        POR = vx[0]
        shipId = vx[1].split(".")[0]
        id = str(POR) + "_" + str(shipId)
        n = vx[2]
        if id in N_R.keys():
              N_R[id].append(n)
for key,value in N_R.items():
    xc = key
    if len(value)==2:
        v1 = value[0]
        v2 =value[1]
        xc=xc+" "+str(v1)+" "+str(v2)+"\n"
        info.write(str(xc))
print(len(N_R))

分港口统计第一梯队的船对应的第二梯队的船舶的ID

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
one_two = open("C:/Users/ouc/Desktop/twoGroupShip/one_twoGroupShips100.txt")
postship = open("C:/Users/ouc/Desktop/firstGroupShip/thepostship100.txt")
page = one_two.readlines();
page2 = postship.readlines()

# 首先将第一港口的船号 与 第二港口的船号 存储到一map里面
oe_map = {}
for x in page:
    sx = x.strip("\n").split(" ")
    if len(sx) >2:
        oe_map[sx[0]] = sx[1:]
    else:
        oe_map[sx[0]] = []
# 读取港口与第一梯队的船对应关系
for y in page2:
    sy = y.strip("\n").split(" ")
    port = sy[0]
    oe_ship = sy[1:]
    tt = open("C:/Users/ouc/Desktop/one_two_100ship/"+str(port)+".txt","a")
    for oo_s in oe_ship:
        if oo_s in oe_map.keys():
            ss = str(oo_s)+" "
            for ui in oe_map[oo_s]:
                ss= ss+" "+str(ui)
            ss=ss+"\n"
        tt.write(ss)

港口与第一梯队的船以及第二梯队的船的对应关系

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
ship = open("/Volumes/17854111053/期末复习/result/re50.txt")
Treejson = open("/Users/gorge/Desktop/range50.json","a")
one_twoID = open("/Volumes/17854111053/期末复习/twoGroupShip/one_twoGroupShips50.txt")
xx = ship.readlines()
count = 0;
port_shipc = {}
for i in range(47):
    port_shipc[str(i + 1)] = 0
ship_portc = {}
for i in range(47):
    ship_portc[str(i + 1)] = []

for x in xx:
    v = x.strip("\n").split(" ")
    if len(v) > 3 and v[2] != '12960':
        port_shipc[v[0]] = (port_shipc[v[0]] + 1)
        ship_portc[v[0]].append(str(v[1]).split(".")[0])
#创建一个第一梯队的船号对应的第二梯队的船的船号
yy = one_twoID.readlines()
one_twoDir = {}
for iy in yy:
    iyy = iy.strip("\n").split(" ")
    one_twoDir[iyy[0]]=[]
    if len(iyy)>1:
        iyy2 = iyy[1:]
        for sd in iyy2:
            one_twoDir[str(iyy[0])].append(str(sd))
# 主要使用的是根据ship_portc进行拼接
# 首先创建一个父类 map
father = {}
father["children"] = []
father["name"] = "港口"
# 先将港口的名字拼接进去
for i in range(47):
    oil = {}
    oil["children"] = []
    oil["name"] = "港口"+str(i+1);
    for sx in ship_portc[str(i + 1)]:
        oil2 = {}
        oil2["children"] = []
        for sds in one_twoDir[str(sx)]:
            oil3 = {}
            oil3["children"] = []
            oil3["name"]=str(sds)
            oil2["children"].append(oil3)
        oil2["name"] = str(sx);
        oil["children"].append(oil2)
    father["children"].append(oil)
#print(father)
Treejson.write(str(father))

分港口统计第二梯队的船的数量

1
2
3
4
5
6
7
8
9
10
11
12
13
nd = {}
dfs = open("C:/Users/ouc/Desktop/twonumber100.txt","a")
for i in range(1,48):
    ss = open("C:/Users/ouc/Desktop/one_two_100ship/one_two_100ship/"+str(i)+".txt")
    page = ss.readlines()
    nd[i]=0;
    for x in page:
        ll = x.strip("\n").split(" ")
        if len(ll)>1:
            llx = ll[1:]
            nd[i]=len(llx)
    dfs.write(str(i)+" "+str(nd[i])+"\n")
print(nd)

聚类算法

一维数据进行聚类，使用设置阈值间隔的方式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import numpy as np
from pandas import Series
def threshold_cluster(Data_set,threshold):
    #统一格式化数据为一维数组
    stand_array=np.asarray(Data_set).ravel('C')
   # print(stand_array)
    #为一维数组赋值上索引
    stand_Data=Series(stand_array).sort_values()
    #print(stand_Data)
    # 创建索引数组 创建分类结果的数组
    index_list,class_k=[],[]
    # 判断函数是不是还有数据，如果有数据的话就可以继续便利了
    while stand_Data.any():
        # 如果只有一个数据的话  就直接将数据拼接到 index_list 以及 class_k中
        if len(stand_Data)==1:
            index_list.append(list(stand_Data.index))
            class_k.append(list(stand_Data))
            # 将数据stand_Data进行清空
            stand_Data=stand_Data.drop(stand_Data.index)
        # 如果多与两个数据的话 就执行下面的代码
        else:
            # 先将第一个数据的index 放置到 第一个分类中
            class_data_index=stand_Data.index[0]
            #print(class_data_index,"---")
            # 将第一个的数据 放置到第一个分类中
            class_data=stand_Data[class_data_index]
            # 将已经分类的数据删除掉
            stand_Data=stand_Data.drop(class_data_index)
            # 判断已经分类的数据与剩下没有分类数据的差值 然后将差值 转化为一维数组
            if (abs(stand_Data-class_data)<=threshold).any():
                args_data=stand_Data[abs(stand_Data-class_data)<=threshold]
                stand_Data=stand_Data.drop(args_data.index)
                index_list.append([class_data_index]+list(args_data.index))
                class_k.append([class_data]+list(args_data))
            else:
                index_list.append([class_data_index])
                class_k.append([class_data])
    return index_list,class_k
file = open("/Users/gorge/Desktop/thepostship100number.txt")
Data_set = []
page = file.readlines()
for sc in page:
    df= sc.strip("\n").split(" ")
    k = df[0]
    v = df[1]
    Data_set.append(int(v))
index_list, class_k = threshold_cluster(Data_set, 5)
print(index_list)
print(class_k)
第二梯队的船相遇了多少第一梯队的船，与第一梯队的船相遇了多少次 熟悉度是怎么样的

  • 第二梯队的船相遇第二梯队的船多少次

    • 首先是祛除重复的统计出第二梯队的船，放在 Twoships.txt文件中，然后读取第一梯队与第二梯队对应关系的one_twoGroupShipsx.txt的文件，如果该文件第二列以后的ID包含现在关注的第二梯队的船的ID，我们就将 该第一梯队的id,添加到对应的第二梯队的ID下面。这样就完成了第一步统计第二梯队的船相遇第二梯队的船多少次。

  • 第二梯队的船与第一梯队的船相遇的次数是怎么样的

    • 根据船与船之间相遇关系的文件rx_relation.txt的文件，只要船A的ID以及船B的ID能够与我们上面统计的第二梯队与第一梯队的船对应上（这里的对应上有两种情况 船A的ID属于第一梯队以及船B的id属于第二梯队 或者 船A的ID属于第二梯队以及船B的id属于第一梯队 ）这样我们就将这相遇次数 添加进去。
    • 在这个频率下对应的熟悉度的大小是怎么样的
      • 读取rx_waitRealation.txt这个文件，计算船舶与船舶之间的熟悉度的大小关系。按照与统计第二梯队的船与第一梯队的船的相遇次数的方式进行匹配熟悉度
        1
        2
        3
        4
        5
        6
        7
        8
        9
        10
        11
        12
        13
        14
        15
        16
        17
        18
        19
        20
        21
        22
        23
        24
        25
        26
        27
        28
        29
        30
        31
        32
        33
        34
        35
        36
        37
        38
        39
        40
        41
        42
        43
        44
        45
        46
        47
        48
        49
        50
        51
        52
        53
        54
        55
        56
        57
        58
        59
        60
        61
        62
        63
        64
        65
        66
        67
        68
        69
        70
        ship = open("C:/Users/ouc/Desktop/twoGroupShip/one_twoGroupShips20.txt")
        ship_id = open("C:/Users/ouc/Desktop/twoGroupShip/Twoships20.txt")
        ## 用来存储第二梯队的船相遇第一梯队的船的数目
        two_one_nunmber = open("C:/Users/ouc/Desktop/twoGroupShip/two_one_nunmber20.txt","a")
        page_ship = ship.readlines();
        page_id = ship_id.readlines();
        twog = {}
        twoc = {}

        for id in page_id:
            id_1 = id.strip("\n")
            twoc[str(id_1)] = 0
            for x in page_ship:
                if len(x) > 0:
                    x_1 = x.strip("\n").split(" ")[0]
                    xx = x.strip("\n").split(" ")[1:]
                    if xx.__contains__(id_1):
                        k1 = str(id_1)+"_"+str(x_1)
                        twog[k1] =[]
                        twoc[str(id_1)] = twoc[str(id_1)]+1
        for k,v in twoc.items():
            two_one_nunmber.write(str(k)+" "+str(v)+"\n")
        # 第二梯队的船与第一梯队的船相遇的次数是怎么样的
        two_oneship = open("C:/Users/ouc/Desktop/result/r20_relation.txt")
        two_onepage = two_oneship.readlines()
        twog_key = twog.keys()
        for vd in two_onepage:
            line = vd.strip("\n").split(" ")
            x1 = line[0].split(".")[0]+"_"+line[1].split(".")[0]
            x2 = line[1].split(".")[0] + "_" + line[0].split(".")[0]
            n = int(line[2])
            if twog_key.__contains__(x1) or twog_key.__contains__(x2):
                if twog_key.__contains__(x1):
                    twog[x1].append(n)
                if twog_key.__contains__(x2):
                    twog[x2].append(n)

        # 计算相应的熟悉度并进行添加
        ship_wait = open("C:/Users/ouc/Desktop/result/r20_waitRelation.txt")
        info = open("C:/Users/ouc/Desktop/twoGroupShip/r20_two_one_concretely.txt","a")
        xx = ship_wait.readlines()
        T = 12957
        count = 0;
        for x in xx:
            v =x.strip("\n").split(" ")
            x1 = v[0].split(".")[0] + "_" + v[1].split(".")[0]
            x2 = v[1].split(".")[0] + "_" + v[0].split(".")[0]
            if twog_key.__contains__(x1) or twog_key.__contains__(x2):
                if len(v)>2:
                    v = v[2:]
                    SUM = 0
                    if v[0]!='12960':
                        count+=1;
                        for vi in v:
                            SUM+=(((int(vi)**2)/2))
                        if SUM ==0.0:
                            SUM = T-3
                        else:
                            sx =T/SUM
                if twog_key.__contains__(x1):
                    twog[x1].append(sx)
                if twog_key.__contains__(x2):
                    twog[x2].append(sx)

        for k,v in twog.items():
            strx = k.split("_")[0]+" "+k.split("_")[1]
            for ix in v:
                strx =strx+" "+str(ix)
            strx =strx+"\n"
            info.write(strx)

树状图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
ship = open("C:/Users/ouc/Desktop/twoGroupShip/r20_two_one_concretely.txt")
Treejson = open("C:/Users/ouc/Desktop/twoGroupShip/r20.json","a")
xx = ship.readlines()

eds = {}
for xd in xx:
    line = xd.strip("\n").split(" ")
    if eds.keys().__contains__(line[0]):
        fir = {}
        fir[line[1]] = []
        fir[line[1]].append(line[2])
        fir[line[1]].append(line[3])
        eds[line[0]].append(fir)
    else:
        eds[line[0]]=[]
father = {}
father["children"] = []
father["name"] = "The Two Group ship"
#
for k,v in eds.items():
    one = {}
    one["children"] = []
    one["name"] = str(k)
    for scx in v:
        for xsf in scx.keys():
            num = {}
            num["children"] = []
            num["name"] = xsf
            cishu = {}
            cishu["children"] = []
            cishu["name"] = scx[xsf][0]
            num["children"].append(cishu)
            cf = {}
            cf["children"] = []
            cf["name"] = scx[xsf][1]
            num["children"].append(cf)
            one["children"].append(num)
    father["children"].append(one)
print(father)
Treejson.write(str(father))

### 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
//统计船与港口的相遇与不相遇的时间节点  
	    public static List<String> getMeetportTime(List<String> firstline ,List<List<String>> second) {
			//定义一个存储上一个时刻的变量
			double pretime = 0.0;
			//定义一个变量用来 存储船的相遇次数
			int count = 0;
			String secshipid = second.get(0).get(1);
			String firshipid = firstline.get(0);
			System.out.println("港口"+firshipid+"相遇"+secshipid+"getMeetCount() start");
			List<String> lStrings = new LinkedList<String>();
			lStrings.add(firshipid);
			lStrings.add(secshipid);
			//初始化一个用来记录两条船相遇时间 节点的链表[1,2]
			List<String> timenode = new LinkedList<String>();
			//定义一个变量用来存储在船相遇的时候，伴随了几个间隔
			int b = 0;
			//定义一个旗帜来表示这段间隔是不是船开始伴随了
			boolean flag = false;
			//定一个变量来判断表示两个船是不是第一次相遇
			boolean firstmeetTime = true;
			//定义一个变量来存储第一次相遇的时刻
			double firtimeet = 0;
			//定义一个变量用来保存 最后一次相遇的时刻
			double endtimeet = 0;	
			if(second.size()>=4) {
			for(int i=0;i<second.size()-1;i++) {
				//获得第二条船的第i时刻的记录
				 List<String> secondline = second.get(i);
				 //记录一下当前时刻 
				 double curtime = Double.parseDouble(secondline.get(0));
				 //获得港口的第i时刻的 x y  
				 double firstx = Double.parseDouble(firstline.get(0));
				 double firsty = Double.parseDouble(firstline.get(1));
				 //获得第二条船的第i时刻的x y
				 double secondx = Double.parseDouble(secondline.get(2));
				 double secondy = Double.parseDouble(secondline.get(3));
				 //计算一下 两个点之间的距离 与 船舶的二倍的通信半径比较一下
				 double x = Math.pow((secondx-firstx), 2);
				 double y = Math.pow((secondy-firsty), 2);
				 double distance = Math.sqrt(x+y);
				 //设置一下船的通信范围
				 double contactRange = 20*1.87;
				 if(distance<=contactRange) {
					 if(firstmeetTime) {
						 firtimeet = curtime;
						 firstmeetTime = false;
					 }
					 endtimeet = curtime+3;
					 if(curtime-pretime>3) {
					 pretime = curtime;
					 count++;
					 flag= true;
					 b++;
					 }else {
						pretime = curtime;
						flag = true;
						b++;
						//System.out.println("b"+b);
					}
				 }else {
					 flag=false;
				 }
				 if(!flag && b!=0) {
					// System.out.println(flag);
					 timenode.add(String.valueOf("["+(pretime-b*3)+","+pretime+"]"));
					 b=0;
				 }
			}
			}
			if(b!=0) {
				//System.out.println("out"+b);
				timenode.add(String.valueOf("["+(pretime-b*3)+","+pretime+"]"));
			}
			//伴随的次数
			 lStrings.add(String.valueOf(count));
			 //每次伴随的时间节点对
			 lStrings.add(String.valueOf(timenode));
			 //完成伴随需要花费的总的时间
			 lStrings.add(String.valueOf(endtimeet-firtimeet));
			System.out.println("船"+firshipid+"相遇"+secshipid+"次数为："+count);
			System.out.println("getMeetCount() end");
			return lStrings;
		}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import matplotlib.pyplot as plt
filepath = open("C:/Users/ouc/Desktop/result/re20meetTimeNode2.txt")
filepath2 = "C:/Users/ouc/Desktop/picture/"
page = filepath.readlines()
ddct = {}
for t in range(24574560,24587517,3):
    ddct[t] = 0
for line in page:
    lines = line.strip("\n").split(" ")
    if lines[2]!='0':
        xd = lines[3:-1]
        p = filepath2+lines[0]+"_"+lines[1]+'.png';
        for v in xd:
            vv = v.strip("[").strip("]").strip(",").split("_")
            x1 = int(vv[0])
            x2 = int(vv[1])
            for ss in range(x1,x2,3):
                ddct[ss] = 1
            print(vv)
        x = []
        y = []
        for key, val in ddct.items():
            x.append(key)
            y.append(val)
        plt.plot(x, y)
        plt.savefig(p,dpi=300)
        plt.close()

计算两种熟悉度的数值 并进行相关统计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
meet =  open("C:/Users/ouc/Desktop/result/re20meet.txt")
wa = open("C:/Users/ouc/Desktop/result/re20.txt")
ress = open("C:/Users/ouc/Desktop/result/resss.txt","a+")
# 初始化一个meet 次数
n_map = {}
c_map={}
meet_ship = meet.readlines()
for line in meet_ship:
    ls = line.strip("\n").split(" ")
    if str(ls[2])!="0":
        k = str(ls[0])+"_" + str(ls[1].split(".")[0])
        n_map[k] = ls[2]
wa_ship = wa.readlines()
wa_map = {}
for line in wa_ship:
    ls = line.strip("\n").split(" ")
    if str(ls[2])!="12960":
        k = str(ls[0]) + "_" + str(ls[1].split(".")[0])
        fg = ls[2:]
        if len(fg)>0 and n_map.__contains__(k):
            wa_map[k] = []
            for ys in fg:
                wa_map[k].append(ys)
###  第一个公式
def gef(ti,n):
    sum = 0;
    for x in ti:
        sum = sum+int(x);
    T = sum/len(ti)
    f = 0;
    for y in ti:
        f = f+(int(y)-T)**2
    f=f/len(ti)
    return n/(f**0.005);
# 第二个公式
def getsd(loc,T):
    fd = 0
    for i in range(1,len(loc)):
        fd =fd+ ((loc[i]-loc[i-1])**2)/2
    return T/fd
#
for key,item in wa_map.items():
     ti = item
     n = int(n_map[key])
     if len(ti)>0:
         z = gef(ti,n)
         if str(n)!="1":
            sddx = str(key).split("_")[0]+" "+str(key).split("_")[1]+" "+str(z)+" "+str(n)+"\n"
            ress.write(sddx)

# import time
# st = "[3432]]"
# s = st.strip("[").strip("]")
# print(s)
# time.sleep(4)
# print(len(s))

港口聚簇分类后对港口进行重新编号

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def getNewPort():
    port_map = {}
    port_map[1] = 1
    port_map[2] = 2
    port_map[3] = 3
    port_map[4] = 4
    port_map[5] = 4
    port_map[6] = 5
    port_map[7] = 6
    port_map[8] = 7
    port_map[9] = 8
    port_map[10] = 9
    port_map[11] = 9
    port_map[12] = 9
    port_map[13] = 10
    for i in range(14,22):
        port_map[i]=11
    port_map[22] = 12
    port_map[32] = 13
    port_map[24] = 14
    port_map[26] = 14
    port_map[35] = 15
    port_map[34] = 16
    port_map[23] = 17
    port_map[33] = 17
    port_map[43] = 18
    port_map[27] = 19
    port_map[28] = 19
    port_map[29] = 19
    port_map[37] = 19
    port_map[30] = 20
    port_map[36] = 21
    port_map[38] = 22
    port_map[46] = 22
    port_map[39] = 23
    port_map[40] = 24
    port_map[31] = 25
    port_map[41] = 26
    port_map[42] = 26
    port_map[43] = 26
    port_map[44] = 26
    return port_map
ff = getNewPort()
print(ff)

问题

原始文件 大小11G，这11G数据是所有的船的数据都混合在一起的，使用Python 进行按照船的编号进行分文件，将同一条船的数据放在一个txt 文件下，分出来的单条船组成的文件的大小，要比原始文件的大小小很多，程序也没有报错！

划分网格 统计船舶在相应海域的密度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import os
rootpath = "C:/Users/ouc/Desktop/data/New201609/Sample1"
targetpath = open("C:/Users/ouc/Desktop/result/201609s.txt","a+")
ls = os.listdir(rootpath)
arr = [[0 for col in range(1000)] for row in range(0,1000)]
def getY(lon1,lat1,lon2,lat2,lat3):
    k = (lat2-lat1)/(lon2-lon1)
    b = lat2-k*lon2
    return (lat3-b)/k
#设置经纬度的范围
minLat = 25
maxLat = 35
minLon = 120
maxLon = 130
Ilon = (maxLon-minLon)/1000
Ilat = (maxLat-minLat)/1000
print(Ilon)

print(len(arr))
for p in ls:
    path = rootpath+"/"+p;
    page = open(path).readlines()
    for line in page:
        l = line.strip("\n").split(" ");
        lon = float(l[4])/600000
        lat = float(l[3])/600000
        if lon>=120 and lon<=130 and lat<=35 and lat>=25:
            x = int((lon - minLon) / Ilon)
            y = int((lat - minLat) / Ilat)
            lony =0
            if lat>=32.62170 and lat<=35:
                lony = getY(121.39759,34.040604,121.898542,32.62170,lat)
            if lat<32.62170 and lat>=32.02698:
                lony = getY(121.898542,32.62170,122.742584,32.02698,lat)
            if lat<32.02698 and lat>=30.68529:
                lony = getY(122.742584,32.02698,123.836783,30.68529,lat)
            if lat<30.68529 and lat>=28.312783:
                lony = getY(123.836783,30.68529,122.613561,28.312783,lat)
            if lat < 28.312783 and lat >= 25.625886:
                lony = getY(122.613561,28.312783, 120,25.625886, lat)
            if lony < lon:
                print(lony,":",lon)
                arr[y][x] = arr[y][x] + 1;

mx = 0
for i in arr:
    li =""
    for j in i:
        li=li+(" "+str(j))
        if j>mx:
            mx = j
    targetpath.write(li+"\n")
print(mx)
统计船所属的港口